knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, 
                      fig.height = 4, fig.width = 7)

library(terra)
library(sf)
library(data.table)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))

Summary

Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.

Methods

Summarize functional entity membership per cell across IUCN and AquaMaps

Included species are:

  • coded for functional entity membership
  • mapped in one (or both) of the species distribution mapsets
  • included the trait-based vulnerability set

Load species ranges and summarize to functional entities per cell

Here we use the IUCN species range maps rasterized to 10 km Mollweide, and AquaMaps species distributions reprojected to 10 km Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥ 10. These map filenames are already set up in the assemble_spp_info_df() function.

Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species for later calculation into functional richness, redundancy, vulnerability. Break the process into chunks to avoid overloading memory.

process_spp_to_fe <- function(fe, spp_maps) { ### fe <- fe_vec[1]
  dt_out <- spp_maps %>%
    filter(fe_id == fe) %>%
    data.table() %>%
    .[ , .(n_spp_fe = length(unique(species))),
       by = .(cell_id, fe_id)]
  return(dt_out)
}
fe_summary_file <- here_anx('func_entities', 'fe_cell_summary.csv')
# unlink(fe_summary_file)
if(!file.exists(fe_summary_file)) {

  spp_info_df <- assemble_spp_info_df(fe_only = TRUE, vuln_only = TRUE)
  ### table(spp_info_df %>% select(species, len) %>% distinct() %>% .$len)
  # huge   lg  med   ml   sm tiny  vlg  vsm 
  #  406 1370 3595 2555 4179 2174  906 5959 

  ### Chunk out and save smaller files to avoid loading in EVERY spp map
  ### all at once!
  n_chunks <- 10
  spp_to_chunks <- spp_info_df %>%
    select(species) %>%
    distinct() %>%
    mutate(chunk = rep(1:n_chunks, length.out = n()))
  
  spp_clean <- spp_info_df %>%
    select(species, fe_id, map_f, taxon) %>%
    distinct() %>%
    left_join(spp_to_chunks, by = 'species')
  
  spp_fes <- spp_clean %>%
    select(species, fe_id) %>%
    distinct()

  tmp_chunk_filestem <- here('tmp/fe_summary_chunk_%s.csv')
  for(i in 1:n_chunks) { ### i <- 1
    tmp_chunk_f <- sprintf(tmp_chunk_filestem, i)
    if(file.exists(tmp_chunk_f)) {
      message('Temp file ', basename(tmp_chunk_f), ' exists... skipping!')
      next()
    }
    message('Assembling species maps and binding FEs for chunk ', i, ' of ', n_chunks, '...')
    spp_chunk <- spp_clean %>% filter(chunk == i)
    chunk_spp_maps <- collect_spp_rangemaps(spp_vec  = spp_chunk$species, 
                                            file_vec = spp_chunk$map_f,
                                            idcol = 'species',
                                            ) %>%
      dt_join(spp_fes, by = 'species', type = 'left') %>%
      distinct()
    
    fe_vec <- chunk_spp_maps$fe_id %>% unique() %>% sort()
    
    message('... summarizing FR map for chunk ', i, ': ', length(fe_vec), ' FEs in ', 
      nrow(chunk_spp_maps), ' cell observations...')
    chunk_fe_sumlist <- parallel::mclapply(fe_vec, mc.cores = 40,
                                           FUN = process_spp_to_fe,
                                           spp_maps = chunk_spp_maps)
    if(check_tryerror(chunk_fe_sumlist)) {
      stop('Try-errors detected in FE summary loop!')
    }
    
    message('... binding chunk and writing chunk summary to ', basename(tmp_chunk_f), '...')
    chunk_fe_summary <- chunk_fe_sumlist %>% data.table::rbindlist()
    fwrite(chunk_fe_summary, tmp_chunk_f)
  }
  
  message('... reading temp chunks, binding, and summarizing by fe and cell...')
  fe_tmp_fs <- list.files(dirname(tmp_chunk_filestem), 
                          pattern = 'fe_summary_chunk', full.names = TRUE)
  fe_summary_list <- parallel::mclapply(fe_tmp_fs, FUN = fread, mc.cores = 10)
  fe_summary <- rbindlist(fe_summary_list) %>%
    .[ , .(n_spp_fe = sum(n_spp_fe)),
       by = .(cell_id, fe_id)]
  
  message('... writing summary to ', basename(fe_summary_file), '...')
  fwrite(fe_summary, fe_summary_file)
  ### unlink(fe_tmp_fs)
}

Calculate functional diversity metrics

Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.

This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.

  • spp_fe_df joined to am_spp_cell;
  • group by cell and functional group;
  • calculate number of spp per functional group;
  • group by cell;
  • calculate functional entity metrics -
    • Species richness
    • Functional entity richness
    • Functional vulnerability (% of FEs represented by a single spp)
    • Weighted functional vulnerability (new: vuln of an FE is \(\frac{1}{2^{n-1}}\), take average over all FE)
    • Functional redundancy (mean spp count per FE)
    • Functional over-redundancy (sum of spp per FE above mean, divided by total number of spp)
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')

### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
  # system.time({
  fe_sum_total <- data.table::fread(fe_summary_file)

  ### there are 6.5 million cells.  Chop into 500k instances and mclapply to summarize
  # cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
  chunk_size <- 500000
  n_chunks <- 6.5e6 / chunk_size

  message('Calculating species richness map...')
  nspp_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
                 FUN = function(n) { ### n <- 1
                   cell_id_min <- (n - 1) * chunk_size + 1
                   cell_id_max <- n * chunk_size
                   nspp_sum <- fe_sum_total %>%
                     filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                     data.table() %>%
                     .[ , .(n_spp = sum(n_spp_fe, na.rm = TRUE),
                            n_fe  = n_distinct(fe_id)),
                        by = 'cell_id']
                   }) 
  
  if(check_tryerror(nspp_list)) stop('Try errors detected!')
  nspp_df <- nspp_list %>%
    data.table::rbindlist()

  message('Calculating functional vulnerability map...')
  fv_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
                 FUN = function(n) { ### n <- 1
                   cell_id_min <- (n - 1) * chunk_size + 1
                   cell_id_max <- n * chunk_size
                   f_vuln <- fe_sum_total %>%
                     filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                     oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
                     data.table() %>%
                     .[ , .(f_vuln  = sum(n_spp_fe == 1) / n_fe, 
                            f_wvuln = mean(0.5^(n_spp_fe - 1)),
                            f_red   = mean(n_spp_fe)),
                     by = .(cell_id, n_fe, n_spp)]
                 }) 
  if(check_tryerror(fv_list)) stop('Try errors detected!')
  f_vuln_df <- fv_list %>%
    data.table::rbindlist()

  message('Calculating functional overredundancy map...')
  f_overred_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
              FUN = function(n) { ### n <- 1
                cell_id_min <- (n - 1) * chunk_size + 1
                cell_id_max <- n * chunk_size
                f_overred <- fe_sum_total %>%
                  filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                  oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
                  data.table() %>%
                  .[ , f_over := ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)] %>%
                  .[ , .(f_overred = sum(f_over) / sum(n_spp_fe)),
                     by = 'cell_id']
              })
  if(check_tryerror(f_overred_list)) stop('Try errors detected!')
  f_overred_df <- f_overred_list %>%
    data.table::rbindlist()

  message('Combining functional diversity metrics maps...')
  fe_metrics_df <- f_vuln_df %>%
    oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full')
  
  write_csv(fe_metrics_df, fe_metrics_file)
}

Map functional diversity metrics

Create and save rasters of functional diversity metrics

fe_metrics_df <- data.table::fread(fe_metrics_file)
### note that while there are 3,684,273 cells in this dataframe,
### 4273 of those are the Caspian Sea, so when those are dropped,
### it results in 3,680,000 exactly... ugh, weird

n_spp_rast   <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast    <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast  <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast   <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')

writeRaster(n_spp_rast,   here('_output/func_entities/n_spp.tif'),   
            overwrite = TRUE)
writeRaster(n_fe_rast,    here('_output/func_entities/n_fe.tif'),    
            overwrite = TRUE)
writeRaster(f_vuln_rast,  here('_output/func_entities/f_vuln.tif'), 
            overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'), 
            overwrite = TRUE)
writeRaster(f_red_rast,   here('_output/func_entities/f_red.tif'),   
            overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'), 
            overwrite = TRUE)

Examine distributions of each raster.

Plot each raster. NOTE: These no longer include spp that are missing from vulnerability calculations - so the maps in _output/nspp_maps (vuln \(\cap\) FE) should be identical to these in output/func_entities.